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SUMMARY 

This  report  describes  the  formulation  of  a  three-dimensional  penalty 
element  that  uses  the  bulk  stress  at  each  node  as  a  degree  of  freedom,  in 
addition  to  the  usual  u,  v,  and  w  displacements.  The  penalty  function  method 
is  used  to  define  a  modified  energy  functional  that  is  constrained  by  the 
equation  for  bulk  stress,  and  a  solution  is  obtained  by  specifying  stationarity 
of  this  modified  functional.',  Edge-cracked  and  centre-cracked  panels  are 
modelled  in  order  to  investigate  the  practical  implementation  of  the  penalty 
element  formulation  and  to  compare  the  results  with  those  obtained  using 
standard  methods.  /;  ^  J —  v,  ,  /  /}  /T**"  —7.  - 
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1.0  INTRODUCTION 


In  recent  years  the  thermoelastic  effect  has  been  greatly  exploited  as  a  non-contact 
method  of  experimental  stress  analysis.  A  device  known  as  SPATE  (Stress  Pat¬ 
tern  Analysis  by  measurement  of  Thermal  Emission)  has  been  used  extensively  at 
the  Aeronautical  Research  Laboratory  to  analyze  many  different  components  to 
gain  data  on  the  surface  stress  distributions.  The  SPATE  apparatus  accomplishes 
its  analysis  by  measuring  the  thermal  emission  of  the  surface  of  the  component 
which  is  being  sinusoidally  loaded  about  some  mean  load.  The  thermal  emission 
is  proportional  to  the  surface  temperature  and,  because  the  thermal  fluctuations 
are  proportional  to  the  changes  in  bulk  stress  ( ozz  +  oyy  +  o2Z),  a  strain  gauge 
can  be  used  to  calibrate  the  output  from  SPATE  in  order  to  determine  the  actual 
bulk  stress. 

Although  the  use  of  SPATE  for  the  analysis  of  components  is  a  powerful  and 
versatile  technique,  the  tensorial  nature  of  stress  is  ignored  since  temperature  and 
bulk  stress  are  scalar  quantities.  This  places  a  significant  limitation  on  the  ability 
of  SPATE  data  to  provide  a  quantitative  measure  of  the  stress  field,  except  in 
circumstances  where  it  is  known  that  the  stress  field  is  primarily  uni-directional. 
Another  disadvantage  is  that  no  direct  information  concerning  the  nature  of  the 
internal  stress  field  of  a  three-dimensional  structure  can  be  obtained  from  the 
simple  surface  analysis  performed  by  SPATE. 

It  is  well  known  that  the  stress  field  existing  in  a  component  must  satisfy 
equilibrium  and  any  imposed  boundary  conditions.  For  two-dimensional  problems, 
Ryall  and  Wong  [l]  have  shown  that  when  additional  stress-related  information  is 
available,  such  as  experimentally  obtained  bulk  stress  data  from  a  SPATE  analysis, 
it  is  no  longer  necessary  to  have  a  complete  knowledge  of  the  boundary  conditions 
governing  the  problem  in  order  to  obtain  a  useful  solution.  Their  method  involved 
the  use  of  an  assumed  stress  potential  $,  with  unknown  coefficients,  from  which 
the  bulk  stress  could  be  calculated  analytically.  The  values  of  the  unknown  co¬ 
efficients  were  determined  by  formulating  the  problem  using  a  best-fit  criterion 
where  the  difference  between  the  analytical  expression  for  the  bulk  stress  and  the 
experimentally  observed  bulk  stresses  was  minimized  in  the  least  squares  sense. 

This  report  describes  the  formulation  of  an  alternative  finite  element  approach 
using  a  three-dimensional  penalty  element  that  uses  the  bulk  stress  at  each  node 
as  a  degree  of  freedom,  in  addition  to  the  usual  u,  v,  and  w  displacements.  The 
penalty  function  method,  as  described  in  [2],  will  be  used  to  define  a  modified 
energy  functional  that  is  constrained  by  the  equation  for  bulk  stress,  a  solution 
being  obtained  by  specifying  stationarity  of  this  functional. 

This  approach  leads  to  an  interesting  possible  application  where  a  finite  el¬ 
ement  model  of  a  component  is  created,  and  the  bulk  stress,  as  obtained  by  an 
experimental  technique  such  as  SPATE,  is  prescribed  at  various  points.  The  re¬ 
sulting  system  can  then  be  solved,  making  use  of  any  other  available  data  concern¬ 
ing  the  applied  loads  and  boundary  conditions,  thus  determining  the  full  three- 
dimensional  stress  field  in  the  component.  In  problems  where  ill-conditioning  is 


encountered,  such  as  adhesively  bonded  repairs,  it  may  also  be  possible  to  use  the 
measured  bulk  stress  field  as  a  means  of  improving  the  conditioning,  thus  driving 
the  solution  to  a  more  reliable  answer. 


2.0  MATHEMATICAL  FORMULATION 

In  the  finite  element  method  the  displacements  u  at  any  point  within  an  element 
may  be  written  as 

u  =  Na,  (2.1) 

where  the  components  of  N  are  prescribed  functions  of  position  and  a  represents 
the  vector  of  nodal  displacements  for  a  particular  element.  The  vector  a  may  be 
written  in  terms  of  the  n  nodal  points  as 

»i 

a2 

a„ 

and  if  each  a,  is  allowed  horizontal,  lateral,  and  vertical  displacement  degrees  of 
freedom  then  we  can  write 
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v i  >  ■  (2-3) 
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With  the  displacements  known  at  all  nodal  points,  then  by  using  the  linear 
theory  of  elasticity  the  strains  can  be  written  as 
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and  the  relationship  between  stress  and  strain  is  of  the  form 


(2.5) 
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(2.6) 


where  D  is  an  elasticity  matrix  containing  the  appropriate  material  properties. 

We  will  now  introduce  the  definitions  of  the  bulk  stress  n b,  the  bulk  strain  e, 
and  the  bulk  modulus  of  the  material  k: 


&b  —  &xz  "f  &yy  O zz 

(2.7) 

e  =  ez  -f  ey  +  ez 

(2.8) 

E 

(2.9) 

K  =  - — 

1  —  2i/ 

where  axz  oyy  azl  are  the  principal  stresses,  and  ex  ey  ez  are  the  principal  strains, 
E  is  Young’s  modulus,  and  i>  is  Poisson’s  ratio. 

The  equation  relating  the  bulk  stress  Oj,  to  the  bulk  strain  e  is 


(2.10) 

(2.11) 

(2.12) 

and  the  bulk  modulus  vector  K  is  defined  as 
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(2.13) 


Equation  (2.12)  may  be  written  in  the  form 

ab  -  Kre  =  0 


and  we  note  that  the  bulk  stress  ab  is  a  scalar  quantity. 

We  now  proceed  to  define  an  energy  functional  n  which,  in  the  absence  of 
body  forces,  is  a  combination  of  the  standard  energy  functional  Ila  due  to  elastic 
deformations,  and  a  penalty  function  term  n„  which  involves  the  bulk  stress  ab, 
such  that 


n(e,ot)  =  na  +  n„ 


n a  =  ^J  erDedV  -  I  utT  dS 
n *  =  1/  (cb-Krey  dv. 


and  T  is  the  surface  traction  vector  {Tx,Tv,Tz}t . 

Now,  by  substituting  for  u,  the  strains  e  can  be  written  in  terms  of  the  nodal 
displacements  as  follows 


e  =  BNa. 


Hence,  the  standard  energy  functional  IIa  may  be  expressed  as 


na  =  2aTRa  ~  aTfr 


where  the  stiffness  matrix  K  and  the  force  vector  fT  are  given  by 


k  =  /n’ 


BtDBN  dV 


and  the  penalty  functional  IT,,,  becomes 


n„  =  f  J  (ob  -  2«7(,/ctc  +  *Tc/cTe)  dV  (2.22) 
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and  P  is  a  penalty  parameter  whose  nominal  value  is  assigned  such  that  the  order 
of  magnitude  of  the  two  functionals  Ila  and  na  is  approximately  the  same.  In  our 
particular  case,  it  will  be  shown  later  that  the  penalty  parameter  can  be  specified 
to  be  equal  to  either  1/E  or  1/k  without  significant  change  in  the  solution. 

Let  us  now  define  <rb  in  terms  of  a  shape  function  M  and  the  values  of  the 
bulk  stress  at  the  nodal  points,  the  latter  denoted  by  the  bulk  stress  vector  b,  such 
that 


ob  =  MTb, 

(2.28) 

and  by  noting  that  MTb  =  bTM  we  obtain  a  new 
functional 

expression  for  the  penalty 

n,r  =  ~  J  bTMMTb  dV  -  P  j  bTM*TBN  dV  a  +  yaT  Kttaa  (2.29) 

=  ^bTKwb  -  PbTK(,aa  +  ^aTKaaa 
«  i 

(2.30) 

where 

KM  =  j  MMt  dV 

(2.31) 

Kba  =  j  M*TBNdF. 

(2.32) 
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Hence  the  complete  expression  for  the  energy  functional  fl  is 


Tl(a,b)  =  ^aTKa  -  aTfT  +  £bTKwb  -  PbTK6aa  +  ^aTKaaa  (2.33) 
z  z  z 


where  n  is  now  a  function  of  the  displacement  degrees  of  freedom  a  and  the  bulk 
stress  degrees  of  freedom  b  specified  at  each  node. 

In  order  to  obtain  a  solution,  we  apply  the  principle  of  stationarity  to  the 
energy  functional  IT,  namely 


da 

(2.34) 

^=0 

db 

(2.35) 

expressions 

I11  =  Ka  -  f_  +  PKaaa  -  PKbab 
d&  1 
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(2.36) 

(2.37) 


Hence,  the  stationarity  condition  leads  to  a  set  of  linear  equations,  which  can  be 
expressed  in  matrix  form  as 


K 
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+  P 


Kaa 

-Kta 


-Kjtt 


)  {b}  =  {  0  }  <2  38> 


We  note  that  the  term  involving  K  is  just  the  standard  finite  element  formulation 
of  the  stiffness  matrix,  while  the  matrix  term  multiplied  by  the  penalty  parameter 
P  originates  from  the  penalty  element  formulation  for  the  bulk  stress  degrees  of 
freedom. 


3.0  VALIDATION  OF  PENALTY  ELEMENT  FORMULATION 

In  order  to  check  the  theoretical  formulation  of  the  penalty  element  approach,  as 
well  as  its  practical  implementation  into  a  finite  element  package,  two  demonstra¬ 
tion  problems  will  be  analyzed.  These  consist  of  three-dimensional  finite  element 
models  of  an  edge-cracked  panel  and  a  centre-cracked  panel,  from  which  it  is  pos¬ 
sible  to  determine  bulk  stresses  and  displacements  at  various  nodes  in  the  finite 
element  mesh.  Figures  1  and  2  show  the  general  dimensions  of  the  edge-cracked 
and  centre-cracked  panels  that  were  studied. 

The  test  problems  were  analyzed  using  the  PAFEC  finite  element  package 
running  on  the  ELXSI  6400  computer  at  ARL.  The  penalty  element  formulation 
was  incorporated  into  the  PAFEC  program  by  re-writing  the  code  that  deals  with 
20-noded  orthotropic  brick  elements  in  order  to  incorporate  the  penalty  element. 
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The  penalty  element  was  assigned  an  element-type  number  of  37150,  and  this 
element  ,jpe  overwrites  the  PA KEC  isotropic  element  of  the  same  name.  The 
ortt'o*  ropic  element  type  was  chosen  for  this  modification  because  the  required 
changes  were  easier  to  code  than  for  an  isotropic  element. 


3.1  Edge-cracked  and  centre-cracked  panels 

Figure  3  shows  the  finite  element  mesh  used  for  modelling  the  edge-cracked  panel, 
and  it  details  the  major  dimensions,  including  those  of  the  simulated  crack.  Due 
to  symmetry  only  1/4  of  the  structure  is  shown,  and  the  full  model  represented  by 
this  mesh  can  be  obtained  by  reflection  in  the  X-Z  and  X-Y  planes.  Therefore, 
the  nodes  along  the  z-axis  correspond  to  interior  points  in  the  plate.  The  loading 
applied  to  the  finite  element  model  consists  of  a  uniform  prescribed  e-displacement 
of  0.01  mm  along  its  top  surface  (see  Figures  1  and  2),  and  it  should  be  noted  that 
restraints  are  applied  on  that  face  to  prevent  any  strains  in  the  i-direction. 

The  centre-cracked  panel  was  modelled  by  applying  additional  restraints  to 
the  finite  element  model  of  the  edge-cracked  panel,  making  use  of  the  existing 
mesh.  Hence,  with  the  appropriate  restraints,  the  finite  element  mesh  shown  in 
Figure  3  also  corresponds  to  1  /8  of  the  full  centre-cracked  panel,  and  the  full  model 
can  be  obtained  by  reflecting  the  mesh  in  the  X-Y ,  X-Z  and  Y-Z  planes.  The 
loading  for  the  centre-cracked  panel  was  identical  to  that  used  for  the  edge-cracked 
panel. 

The  material  chosen  for  the  two  different  types  of  panel  corresponds  to  an 
epoxy  adhesive,  whose  material  properties  are  presented  in  Table  1.  Although 
these  properties  are  input  as  for  an  orthotropic  material,  they  are  chosen  so  as 
to  make  the  material  behave  in  an  isotropic  manner.  The  data  are  given  for  two 
different  values  of  Poisson’s  ratio,  where  i/  =  0.30  corresponds  to  the  values  used 
for  the  centre-cracked  panel,  and  1/  =  0.35  corresponds  to  the  edge-cracked  panel. 
Initially  some  of  the  finite  element  calculations  were  carried  out  using  u  =  0.35, 
but  in  order  to  match  the  Poisson’s  ratio  used  for  obtaining  the  stress  intensity 
factor  for  a  centre-cracked  panel  [3],  it  was  also  necessary  to  use  a  Poisson’s  ratio 
of  u  =  0.30. 

Figure  4  shows  a  typical  deformed  finite  element  mesh  that  was  obtained 
for  the  edge-cracked  panel  when  it  was  loaded  by  a  uniform  displacement.  The 
deformation  of  the  mesh  in  the  vicinity  of  the  crack  tip  can  be  seen  clearly. 


3.2  Bulk  stress  calculations 

The  edge-cracked  and  centre-cracked  panels  were  analyzed  using  both  the  standard 
and  penalty  element  formulations.  For  the  standard  element  case,  bulk  stresses 
were  calculated  at  a  number  of  nodes  ahead  of  the  crack  tip  by  simply  summing 
the  principal  stresses  according  to  Equation  (2.7).  Since  the  nodes  chosen  for  this 
comparison  were  placed  at  element  boundaries,  the  bulk  stresses  were  calculated 
as  the  average  of  the  results  obtained  for  adjacent  elements. 
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Tables  2  and  3  present  the  bulk  stresses  calculated  for  the  edge-cracked  and 
centre-cracked  panels.  When  using  the  penalty  element  formulation,  two  different 
values  of  the  penalty  parameter  were  tried,  the  first  corresponding  to  P  =  l/E 
and  the  second  to  P  =  1/it.  The  data  presented  in  these  two  tables  indicates 
that  there  is  very  good  agreement  between  the  standard  and  the  penalty  element 
formulations.  In  areas  where  there  is  a  rapid  stress  variation,  the  answers  generally 
agree  to  within  7%,  while  the  average  difference  is  about  4%.  Note  that  the  results 
obtained  for  the  two  different  values  of  the  penalty  parameter  are  also  in  very  good 
agreement,  indicating  that  the  solution  is  not  very  sensitive  to  changes  in  the  value 
of  the  penalty  parameter  when  it  is  of  the  order  of  magnitude  0(1/E).  Also  note 
that  there  is  a  stress  singularity  at  node  12,  since  this  is  right  at  the  crack  tip, 
and  consequently  the  finite  element  results  are  inaccurate  at  that  location;  the 
results  for  node  12  were  included  only  to  enable  a  qualitative  comparison  to  be 
made  between  the  standard  and  penalty  element  formulations. 


3.3  Stress  intensity  factor  calculations 

In  order  to  give  the  required  \fr  displacement  behaviour  for  near  crack  tip  points 
in  the  vicinity  of  the  crack  front,  special  elements  were  used  that  had  their  midside 
nodes  shifted  to  the  quarter  point  position.  The  value  of  the  stress  intensity  factor 
Kj  near  the  crack  front  can  be  determined  from  the  opening  displacement  of  the 
crack  near  the  tip  using  the  plane  strain  equation 


SE  j  U 

"4(1-*',>V'(1-S) 


(3.1) 


where 

6  is  the  displacement  in  the  y-direction  of  a  point  lying  in  the  X-Z  plane 
near  the  crack  front 

r  is  the  perpendicular  distance  of  the  point  from  the  crack  front 

l  is  the  crack  length  for  an  edge-cracked  panel,  and  is  the  crack  half-length 
for  a  centre-cracked  panel 

E  is  Young’s  modulus 

v  is  Poisson’s  ratio. 

Table  4  presents  values  of  the  stress  intensity  factor  Ki  as  obtained  using  both 
the  standard  and  penalty  element  formulations,  and  compares  them  to  the  values 
of  K[  given  by  Rooke  it  Cartwright  [3]  for  both  an  edge-cracked  and  centre-cracked 
panel.  The  calculations  were  performed  for  an  applied  v-displacement  of  0.01  mm, 
and  it  is  seen  that  the  two  different  finite  element  formulations  give  results  that  are 
in  good  agreement  with  those  obtained  from  [3],  Note  that  the  results  obtained 
using  standard  three-dimensional  elements  are  closer  to  those  given  in  [3],  but  in 
any  case  the  maximum  difference  is  only  about  4%,  which  is  not  very  significant 
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when  the  accuracy  of  the  results  from  [3]  are  quoted  as  being  within  an  error 
bound  of  a  similar  magnitude.  However,  the  difference  in  values  of  Ki  calculated 
using  standard  and  penalty  elements  is  only  about  1.5%,  which  is  very  small. 

Tables  5  and  6  present  values  of  crack  opening  displacement  6,  stress  intensity 
factor  Kj,  and  reaction  force  Fy  obtained  for  the  edge-cracked  and  centre-cracked 
panels.  In  order  to  study  the  sensitivity  of  the  penalty  element  solution  to  dif¬ 
ferent  values  of  the  penalty  parameter  P,  a  wide  range  of  values  were  used,  with 
P  =  100  being  the  largest  value  that  could  be  used  without  numerical  problems. 
The  maximum  variation  in  crack  opening  displacement  and  stress  intensity  fac¬ 
tor  experienced  for  the  edge-cracked  panel  was  1.5%,  while  for  the  centre-cracked 
panel  it  was  2.0%. 


4.0  CONCLUSION 

By  using  a  penalty  element  formulation,  the  standard  finite  element  method  has 
been  modified  to  allow  the  calculation  of  bulk  stresses  as  an  additional  degree 
of  freedom  in  the  finite  element  model  of  a  structure.  The  bulk  stress  penalty 
elements  have  been  incorporated  into  the  PAFEC  finite  element  package. 

In  order  to  check  the  practical  implementation  of  the  penalty  element  formu¬ 
lation,  edge-cracked  and  centre-cracked  panels  were  modelled,  and  the  solutions 
obtained  were  compared  with  those  calculated  using  the  standard  finite  element 
method.  Also,  both  the  standard  and  penalty  element  formulations  have  been 
used  to  obtain  crack  opening  displacements,  from  which  the  values  of  Kj  were 
calculated  and  compared  with  those  from  [3]. 

The  results  of  these  comparisons  indicate  that  bulk  stresses  can  be  accurately 
predicted  using  the  penalty  element  method,  and  that  there  is  no  significant  effect 
on  the  quality  of  the  results  obtained  for  the  other  degrees  of  freedom.  It  has  also 
been  found  that  the  present  application  of  the  penalty  element  formulation  is  very 
insensitive  to  the  choice  of  penalty  parameter  P,  and  a  suitable  value  for  P  can 
be  defined  by  using  either  Young’s  modulus  or  the  bulk  modulus,  the  appropriate 
values  being  given  by  P  =  1/E  or  P  =  1/k. 

It  is  envisaged  that  this  methodology  will  subsequently  be  used,  in  conjunc¬ 
tion  with  experimentally  obtained  values  for  the  bulk  stress,  to  develop  a  hybrid 
experimental/numerical  procedure  for  determining  individual  stress  components 
based  on  the  work  presented  in  [1], 
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TABLE  1 


Mechanical  Properties  of  Adhesive  For 
Two  Values  of  Poisson’s  Ratio 


MECHANICAL 

PROPERTY 

(MPa"1) 

CENTRE-CRACKED 

PANEL 

v  =  0.30 

EDGE-CRACKED 

PANEL 

v  ~  0.35 

®xx 

5.2910  x  10-4 

5.2910  x  10-4 

S  Y  Y 

5.2910  x  lO"4 

5.2910  x  10'4 

Szz 

5.2910  x  10"4 

5.2910  x  10“4 

SXY 

-1.5873  x  10~4 

-1.8519  X  10-4 

SYZ 

-1.5873  X  10"4 

-1.8519  X  lO"4 

®2X 

-1.5873  x  10~4 

-1.8519  x  10-4 

shxy 

1.3757  x  10"3 

1.4286  X  10  3 

shyz 

1.3757  x  10"3 

1.4286  x  10“3 

SHZX 

1.3757  x  10-3 

1.4286  X  lO-3 
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TABLE  2 


Bulk  Stress  at  Different  Nodes  Calculated  Using 
Standard  and  Penalty  Elements  with  Different 
Values  of  Penalty  Parameter  P 


EDGE-CRACKED  PANEL 

DISTANCE 

BULK  STRESS 

NODE 

FROM 

(MPa) 

NO. 

CRACK  FACE 

Standard 

Penalty  Element 

(mm) 

Element 

P  =  1  /E 

Basil 

mm 

0.0000 

26.31 

27.28 

28.09 

0.1667 

12.29 

11.47 

11.48 

14 

0.5000 

5.19 

5.11 

5.08 

15 

1.6667 

2.65 

2.75 

2.74 

16 

5.0000 

2.02 

2.03 

2.01 

TABLE  3 

Bulk  Stress  at  Different  Nodes  Calculated  Using 
Standard  and  Penalty  Elements  with  Different 
Values  of  Penalty  Parameter  P 


CENTRE-CRACKED  PANEL 

DISTANCE 

BULK  STRESS 

NODE 

FROM 

(MPa) 

_ 

NO. 

CRACK  FACE 

Standard 

Penalty  Element 

(mm) 

Element 

saozi 

P=  l/«e 

MM 

0.0000 

21.70 

22.72 

23.09 

0.1667 

9.98 

9.48 

9.43 

14 

0.5000 

4.65 

4.38 

4.37 

15 

1.6667 

2.48 

2.57 

2.55 

16 

5.0000 

2.00 

2.02 

2.01 

14 


TABLE  4 


Values  of  Stress  Intensity  Factor  K /  as  Obtained  Using  Different  Methods 


SOLUTION 

METHOD 

STRESS  INTENSITY  FACTOR  Ej 
(MPav/m) 

EDGE-CRACKED 

PANEL 

CENTRE-CRACKED 

PANEL 

Rooke  ic  Cartwright  [3] 

■■ 

0.140 

Standard  Elements 

0.137 

Penalty  Elements  (P  =  1/E) 

0.149 

0.135 

TABLE  5 


Finite  Element  Results  For  Crack  Opening,  Stress  Intensity  Factor,  and 
Reaction  Force  Obtained  With  Different  Values  of  Penalty  Parameter  P 


EDGE-CRACKED  PANEL 

PENALTY 

CRACK  OPENING  6 

STRESS  INTENSITY 

REACTION  Fy 

PARAMETER 

AT  NODE  11 

FACTOR  if/ 

AT  NODE  3 

P 

(m) 

(MPa^/in) 

(N) 

0.000x10+°* 

1.4062  xlO- 6 

0.1508 

■PH 

l.OOOxlO-6 

1.4061  xlO-6 

0.1508 

l.OOOxlO-5 

1.4051  xlO-6 

0.1507 

l.OOOxlO-4 

1.3984  xlO-6 

0.1500 

1.587xl0-4| 

1.3956X10-6 

0.1497 

5.291  xlO-4J 

1.3874X  10-6 

0.1488 

18.314 

l.OOOxlO-2 

1.3778  xlO-6 

0.1478 

18.338 

1.000x10+° 

1.3858X10-6 

0.1486 

18.362 

1.000  xl0+‘ 

1.3909X10-6 

0.1492 

18.413 

1.000X  10+2 

1.4280X10-6 

0.1532 

18.827 

*  standard  finite  element  formulation 
f  corresponds  to  P  =  1/k, 
t  corresponds  to  P  =  1/E 


TABLE  6 


Finite  Element  Results  For  Crack  Opening,  Stress  Intensity  Factor,  and 
Reaction  Force  Obtained  With  Different  Values  of  Penalty  Parameter  P 


CENTRE-CRACKED  PANEL 

PENALTY 

CRACK  OPENING  6 

STRESS  INTENSITY 

REACTION  Fy 

PARAMETER 

AT  NODE  11 

FACTOR  Kj 

AT  NODE  3 

P 

(m) 

(MPay/m) 

(N) 

0.000x10+°* 

1.3211  xlO-6 

0.1367 

MM 

l.OOOxlO-6 

1.3211  xlO-6 

0.1367 

l.OOOxlO-5 

1.3205  xlO-6 

0.1366 

l.OOOxlO-4 

1.3160X10-6 

0.1361 

2.116X  10-4f 

1.3124  xlO-6 

0.1358 

5.291  Xl0-4t 

1.3067  xlO-6 

0.1352 

l.OOOxlO-2 

1.2959  xlO-6 

0.1340 

1.000x10+° 

1.2986  xlO-6 

0.1343 

18.602 

1.000X10+1 

1.2976  xlO-6 

0.1342 

18.604 

1.000X10+2 

1.2966X10-6 

0.1341 

18.542 

*  standard  finite  element  formulation 
t  corresponds  to  P  =  1/k 
I  corresponds  to  P  =  l/E 


16 


a  =  1.6667  mm 
b  =  10.0  mm 
h  =  10.0  mm 


4 


Fig.  1 :  Edge-cracked  panel  showing  main  dimensions. 


Uniform  displacement, 
no  rotation,  u  =  0 


h 


h 


A 


^  T  * 


Y 

A 


2b 


~x  x  ^  ^  x  ; 


a  =1.6667  mm 
b  =  10.0  mm 
h=10.0  mm 


Fig.  2  :  Centre-cracked  panel  showing  main  dimensions. 


Fig.  3:  Finite  element  mesh  used  to  model  the  edge-cracked  and  centre-cracked 
panels. 


o 


X 


Fig.  4  :  Typical  deformed  finite  element  mesh  for  the  edge-cracked  panel  loaded  by 
a  uniform  displacement. 
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